%% production_table1.m - PRODUCTION run for Table 1
%
% This script runs the FULL Monte Carlo simulation for Table 1 of the paper
% with ALL 100,000 replications.
%
% MANUSCRIPT SETTINGS (Table 1):
%   - Sample sizes: T ∈ {20, 40, 80, 200, 400, 800}
%   - ARMA(1,1) errors: ρ₁ = θ₁ ∈ {0, 0.3, 0.6}
%   - Intervention timing: ξ ∈ {0.75, 0.25}
%   - Innovation distribution: df ∈ {5, 15, Inf} (t-dist or Gaussian)
%   - Replications: 100,000
%   - Bandwidth: S_T = ⌊0.75 * T^(1/3)⌋
%   - Bootstrap: NOT NEEDED (saves significant computational time!)
%
% PRODUCTION SETTINGS (this script):
%   - Replications: 100,000 (FULL)
%   - All settings match the manuscript exactly
%
% OUTPUT LOCATION:
%   Results/production/table1/df{5,15,Inf}/
%
% ⚠️  WARNING: This will take several hours to complete!
%     Run test_table1.m first to verify setup.
%
% ESTIMATED RUNTIME:
%   - Without bootstrap: ~6-18 hours (depending on hardware, 3x for df loop)
%   - With bootstrap (if enabled): ~60-150+ hours

clear all
close all
clc

%% Configuration
% Sample sizes and correlation values (exactly as in manuscript)
time_span = [20, 40, 80, 200, 400, 800];
correlation = [0, 0.3, 0.6];  % Manuscript Table 1: [0, 0.3, 0.6]

% Degrees of freedom for innovation distribution
% df = Inf: Gaussian (benchmark), df = 15: moderate tails, df = 5: heavy tails
df_values = [5, 15, Inf];

% Production settings
num_replications = 100000;     % PRODUCTION: full 100,000 replications
run_bootstrap = false;         % Table 1 does NOT need bootstrap (CRITICAL for speed!)

% Output folder (base)
output_base = fullfile(pwd, '..', 'Results', 'production', 'table1');

% Pre-sample types: 'short' (25%), 'long' (75%)
% Note: Table 1 uses ξ ∈ {0.75, 0.25} which corresponds to 'long' and 'short'
pre_sample_types = {'short', 'long'};

%% Display configuration and confirmation
fprintf('========================================\n');
fprintf('PRODUCTION RUN: Table 1 (JAE Manuscript)\n');
fprintf('========================================\n');
fprintf('⚠️  WARNING: FULL PRODUCTION RUN\n');
fprintf('Replications: %d\n', num_replications);
fprintf('Bootstrap: %s (saves hours of computation!)\n', mat2str(run_bootstrap));
fprintf('Sample sizes: %s\n', mat2str(time_span));
fprintf('Correlations (ρ₁=θ₁): %s\n', mat2str(correlation));
fprintf('Degrees of freedom: %s\n', mat2str(df_values));
fprintf('Pre-sample types: %s\n', strjoin(pre_sample_types, ', '));
fprintf('Output base: %s\n', output_base);
fprintf('========================================\n');
fprintf('Total simulation runs: %d\n', ...
    length(df_values) * length(pre_sample_types) * length(correlation) * length(time_span));
fprintf('Estimated time: 6-18 hours\n');
fprintf('========================================\n\n');

% Confirmation prompt
% response = input('Continue with production run? (y/n): ', 's');
% if ~strcmpi(response, 'y')
%     fprintf('Production run cancelled.\n');
%     return;
% end

%% Run simulations
fprintf('\nStarting production run...\n\n');
tic
total_runs = length(df_values) * length(pre_sample_types) * length(correlation) * length(time_span);
current_run = 0;

for df = df_values
    % Create df-specific output folder
    if isinf(df)
        df_folder = fullfile(output_base, 'dfInf');
        df_str = 'Gaussian';
    else
        df_folder = fullfile(output_base, sprintf('df%d', df));
        df_str = sprintf('t(%d)', df);
    end

    if ~exist(df_folder, 'dir')
        mkdir(df_folder);
    end

    fprintf('\n=== Innovation distribution: %s ===\n\n', df_str);

    for pre_sample_type = pre_sample_types
        % Create subfolder for each pre-sample type
        subfolder = fullfile(df_folder, pre_sample_type{1});
        if ~exist(subfolder, 'dir')
            mkdir(subfolder);
        end

        for rho1 = correlation
            for t = time_span
                current_run = current_run + 1;
                run_start = tic;

                % Progress indicator
                fprintf('[%d/%d] Running: df=%s, T=%d, ρ=%g, ξ=%s\n', ...
                    current_run, total_runs, df_str, t, rho1, pre_sample_type{1});

                output_file_name = fullfile(subfolder, ...
                    sprintf('_%g_%d.mat', rho1, t));

                % Call simulation function with optional parameters
                run_simulation_func(t, rho1, output_file_name, pre_sample_type{1}, ...
                    'num_replications', num_replications, ...
                    'run_bootstrap', run_bootstrap, ...
                    'df', df);

                run_time = toc(run_start);
                remaining_runs = total_runs - current_run;
                est_remaining_time = run_time * remaining_runs;

                fprintf('        Completed in %.2f minutes\n', run_time/60);
                fprintf('        Estimated remaining time: %.2f hours\n\n', ...
                    est_remaining_time/3600);
            end
        end
    end
end

elapsed_time = toc;
fprintf('\n========================================\n');
fprintf('PRODUCTION simulation completed!\n');
fprintf('Total time: %.2f hours\n', elapsed_time/3600);
fprintf('Average time per run: %.2f minutes\n', elapsed_time/60/total_runs);
fprintf('========================================\n');
fprintf('Output saved to: %s\n', output_folder_name);
fprintf('========================================\n');
